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SUMMARY 


The  objective  of  this  program  is  the  analysis  of  the  fundamen- 
tal limitations  placed  upon  image  restoration  by  the  presence  of 
noise  and  unavoidable  atmospheric  and  optical  system  degradations. 
Our  approach  is  to  treat  individually  and  in  concert  the  several 
factors  which  limit  image  formation  (direct  problem)  and  further 
limit  image  restoration  (inverse  problem),  to  understand  their 
interrelation  and  order  of  importance.  Within  this  overall  frame- 
work, we  have  treated  three  problem  areas  comprising  the  chapters 
of  this  report:  sums  of  random  amplitudes  individually  governed 

by  a lognormal  probability  density  function,  realizability  con- 
straints on  wavefront  covariance,  and  (partial)  moment  behavior 
of  the  transfer  function  random  process,  including  temporal  ef- 


fects and  those  due  to  random  amplitude  errors  acting  in  conjunc- 
tion with  random  phase  errors. 

The  lognormal  probability  density  permeates  much  of  the  cur- 
rent literature  on  optical  propagation  through  the  turbulent  at- 
mosphere and  has  been  the  source  of  a good  deal  of  controversy  in 
various  contexts  of  the  general  problem.  Sums  of  lognormal  ran- 
dom variables  arise  through  aperture  averaging  and  other  related 
processes.  Because  the  lognormal  density  assumes  central  impor- 
tance in  many  atmospheric  propagation  problems,  certain  aspects 
of  the  behavior  of  (sums  of)  lognormally  distributed  variables 
are  considered  in  Chapter  I.  Particular  results  include  the  der- 
ivation of  an  expression  for  the  characteristic  function  of  the 
lognormal  probability  density,  which  is  then  used  to  calculate 
densities  for  sums  of  lognormally  distributed  random  variables; 
a simplified  proof  that  the  lognormal  probability  density  is  not 
determined  by  its  moments  even  though  they  exist;  an  investigation 
of  the  observed  "permanence"  of  the  lognormal  density  in  terms  of 


I 
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its  asymptotic  behavior,  leading  to  a soundly  based  explanation 
of  the  phenomenon. 

Realizability  conditions  on  the  covariance  of  the  wavefront 
aberration  function  are  briefly  stated  in  Chapter  II;  the  class 
of  admissible  covariance  functions  is  restricted  to  those  repre- 
senting (isotropic)  random  surfaces.  The  commonly  used  Kolomogorov 
five-thirds  spectrum  is  shown  to  violate  the  stated  realizability 
conditions,  but  no  severe  practical  consequence  normally  arises 
since  the  sampling  interval  is  finite. 

Investigations  of  the  statistics  of  the  transfer  function  ran- 
dom process,  predominantly  restricted  in  our  prior  work  to  errors 
of  phase,  have  been  extended  to  include  the  statistical  dependence 
of  the  transfer  function  process  on  time  (and  concurrently  on  in- 
tegration time),  as  well  as  on  jointly  present  random  amplitude 
and  phase  errors.  Results  show  that  the  expected  value  of  the 
time  dependent  transfer  function  does  not  depend  on  integration 
time,  but  that  the  general  higher  moment  behavior  is  collection 
time  dependent.  Random  amplitude  effects  are  considered  within 
the  usual  lognormal  and  normal  assumptions  regarding  the  respec- 
tive amplitude  and  phase  probability  densities,  but  without  the 
added  restriction  of  Independence  between  the  amplitude  and  phase 
random  processes.  In  terms  of  the  expected  transfer  function,  a 
product  of  four  terms  arises:  one  standard  deterministic  part, 

one  part  due  to  phase  errors  acting  alone,  one  part  due  solely  to 
amplitude  perturbations,  and  a cross  term  which  is  of  unit  modulus 
and  depends  for  its  existence  on  the  odd  part  of  the  cross  cor- 
relation between  amplitude  and  phase. 
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ABSTRACT 


This  report  is  divided  into  three  chapters  which  treat  three 
problem  areas  pertinent  to  optical  propagation  through  the  tur- 
bulent atmosphere.  In  the  first  chapter,  several  aspects  of  the 
lognormal  probability  density,  which  permeates  much  of  the  cur- 
rent literature  and  which  has  been  the  soui  -*e  of  a good  deal  of 
controversy  in  various  relevant  contexts,  are  treated.  Particular 
developments  include  the  derivation  and  use  of  an  expression  for 
the  characteristic  function  of  the  lognormal  probability  density, 
a simplified  proof  that  the  lognormal  probability  density  func- 
tion is  not  determined  by  its  moments,  and  an  investigation  of 
the  observed  permanence  of  the  lognormal  distribution  via  an  adap- 
tation of  the  Cramer  asymptotic  formula  for  sums  of  independent 
random  variables  having  nonzero  means.  The  second  chapter  of  the 
report  concerns  realizability  constraints  on  the  covariance  of 
the  wavefront  aberration  function,  which  a~e  briefly  stated;  it 
is  shown  that  the  class  of  admissible  covariance  functions  is 
restricted  to  those  that  represent  a (isotropic)  random  surface. 
The  commonly  used  Kolomogorov  five-thirds  spectrum  is  seen  to 
violate  the  stated  realizability  conditions,  but  the  practical 
significance  of  the  violation  is  not  excessive  due  to  the  limited 
sampling  interval.  In  chapter  III,  the  effects  of  integration 
time-dependence  and  random  amplitude  on  the  optical  transfer  func- 
tion are  treated  In  terms  of  the  moments  of  the  transfer  function 
random  process.  It  is  shown  that  the  expected  value  of  the  time- 
dependent  transfer  function  is  not  a function  of  the  Integration 
time  Interval,  but  that  general  higher  moments  (explicitly,  second 
order  averaged  products)  do  depend  on  collection  time.  Random 
amplitude  effects  are  considered  for  the  case  of  lognormal  ampli- 
tude and  Gaussian  phase  errors,  without  assuming  Independence 
tween  the  amplitude  and  phase  perturbations.  In  terms  of  the  ex- 
pected transfer  function,  a product  of  four  terms  arises:  one 

deterministic  part,  one  part  due  to  phase  errors  acting  alone,  a 
similar  part  due  only  to  amplitude  errors,  and  a cross  term  which 
is  of  unit  modulus  and  which  depends  for  its  existence  on  the  odd 
part  of  the  cross  correlation  between  log  amplitude  and  phase. 
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CHAPTER  I 

SUMS  OF  INDEPENDENT  L06N0RMALLY  DISTRIBUTED 
RANDOM  VARIABLES 
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1.  INTRODUCTION 


The  lognormal  probability  distribution  permeates  much  of  the 
current  literature1  on  optical  propagation  through  the  turbulent 
atmosphere  and  has  been  the  source  of  a good  deal  of  controversy 
in  various  contexts  of  the  general  problem.  The  purpose  of  the 
present  paper  is  to  report  on  three  aspects  of  the  lognormal  prob- 
ability density  function  which  are  pertinent  to  the  atmospheric 
propagation  problem. 

The  three  topics  are: 

1.  The  development  of  an  expression  for  the  characteristic 
function  of  the  lognormal  probability  density  function  and 
the  use  of  the  characteristic  function  to  calculate  sums  of 
lognormally  distributed  random  variables. 

2.  Simplified  proof,  via  the  Carleman  criterion,  that  the  log- 
normal probability  density  function  is  not  determined  by  its 
moments  even  though  they  exist.  A discussion  of  the  conse- 
quence of  this  result  as  regards  theoretical  approaches  to 
optical  wave  propagation  in  random  atmospheres. 

3.  The  "permanence"  of  the  lognormal  is  investigated  via  an 
adaption  of  the  Cramer  asymptotic  formula  for  sums  of  inde- 
pendent random  variables  having  nonzero  means. 
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2.  CHARACTERISTIC  FUNCTION  OF  LOGNORMAL  PDF 

Consider  the  sum 

N 

x - 1 xk  , (1) 

k=l  K 

where  the  are  independent  random  variables  all  having  the 
identical  lognormal  probability  density  function  (PDF) 


f'x  ) = — - — 

K x,  o/5t 


-(log  xk-C)72o 


xk  10 


where 

K = E(log  xk)  = <log  xk> 


(2) 


o2  = var ( log  xk)  = E[(log  xk~^)2]  . (3) 

As  a matter  of  typographic  convenience,  we  will  use  E(*)  and 
<•>  interchangeably.  It  is  somewhat  more  convenient  to  rewrite 
Eq.  2 in  the  form 


where 


-Clog(x/x  )]z/2oz 

f(x)  = — e 
xo/ 2rr 


(4) 


(5) 
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The  probability  distribution  of  the  sum  x,  call  it  f (x),  is 
most  easily  obtained  as  the  Fourier  transform  of  the  N th  power 
of  the  characteristic  function.  The  characteristic  function  of 
f(xk)  is 


itx, 


= 4>  (t) 


■f 


itx. 


f(xk)e 


Ldx, 


so  that 


(6) 


rN(x) 


— 00 


<t>  (t) 

Xk 


(7) 


The  usefulness  of  this  approach  is  flawed  by  the  fact  that  the 
characteristic  function  of  the  lognormal  PDF  is  unknown  (at  least, 
it  has  never  appeared  in  the  literature).  Therefore  we  have 
carried  out  an  analysis  which  yields  the  characteristic  function, 
an  outline  of  the  analysis  is  presented  in  the  next  few  paragraphs. 


Upon  setting  xk/xfl  = y in  Eq.  2,  we  can  write  Eq.  6 in  the 

form 


<f>  (t) 

xk 


.00 


a/ 2tt 


-yz/2 a2 


ix„ tey 
e dy 


(8) 


Note  that  E,  enters  only  through  the  product  t = xQt.  This  means 

that  for  a given  a,  <j>  (t)  need  only  be  evaluated  once  since  dif- 

<y 

ferent  values  of  E,  or  xfl  merely  rescale  the  characteristic  function. 
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The  main  contribution  to  the  integral  comes  from  the  integrand 
evaluated  in  the  vicinity  of  y = 0.  This  suggests  that  we  expand 

ei«y  . eiT  elty  [eityV2  ^76...-, 

= eiT  eiTy  l b (ix)yn  , (9) 

n=0  n 

and  integrate  termwise.  However,  the  explicit  derivation  of  the 
b-coefficients  rapidly  becomes  unmanageable. 


A more  systematic  approach  to  obtaining  the  equivalent  of 
the  b-coefficients  is  possible  provided  we  rewrite  Eq.  6 in  the 
form 


<t>  (t)  = 


ix 
_e 

a/2? 


e-y 

J—  00 


-y2/2a2  ix(e-y-l) 


dy 


(10) 


Now  expand 

exp[ix(ey-y-l)  ] = l pjy  a (ix)yn  , (11) 

n=0 


where 


,n 


in(ix)  = — n exp[elT(ey-y-l)]  Q 
dy  y 


(12) 
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These  derivatives,  although  extremely  tedious  if  attempted  by 

hand,  are  easily  obtained  on  a computer  using  symbol  manipulation 

programs.  We  employed  the  symbol  manipulation  program  available 

at  MIT  (Project  MAC)  and  the  results  for  the  first  twelve  a are 

n 

listed  in  Table  1 where  s = ix. 


When  the  series  in  Eq.  11  is  substituted  in  Eq.  10,  we  en- 
counter the  integrals: 


y2m+l  e-y2/.2a’sinTydy  = (.1}m  ff2»+2  .-tV/2  h (ffT) 


(13) 


2m  -y2/2o2  , 

y e cosiydy 


_ / ,vm  2m+l  -T2a2/2  . , . 

- (-1)  (p)  a e h^  (at)  , 


(1^) 


where 


hn(x)  = Hen(x) 


(-l)n 


xz/2  dn 


-xz/2 


dx 


n 


(i5) 


Note  that  these  Hermite  polynomials  are  different  than  those 
customarily  employed  in  theoretical  physics.  The  first  few  are: 


■ 

! 


1 

s 
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h0(x)  = 1 , 


h 3 ( x)  = x3-3x 
hj(x)  = x , h4(x)  = x"-6x2+3 

h2(x)  - x 2— 1 , hs(x)  = x5-10x3+15x  (16) 

The  final  expression  for  the  characteristic  function  is 


<f>  (t) 

xk 


= eiT  e“T 2° 2/2  j 
n=0 


(jg) 

n! 


(if) 


hn(ar) 


(17) 


It  is  a complex-valued  function  of  t since  f(x^)  is  nonzero  only 
over  the  semi-infinite  interval  (0,°°).  The  series  converges 
very  slowly  (if  at  all)  for  large  a,  but  practically  the  only 
values  of  a encountered  In  random  wave  propagation  are  for  o < 1. 
Although  it  is  not  the  purpose  of  this  paper  to  present  elaborate 
numerical  calculations,  we  ran  a series  of  tests  for  a - .25  vary- 
ing the  number  of  terms  M in  the  series  and  found  that  M = 9 was 

sufficient  to  guarantee  convergence.  A plot  of  <f>  (t)  is  shown 

xk 

in  Fig.  1;  note  the  oscillatory  behavior. 


Given  this  information,  the  numerical  evaluation  of  fN(x)  via 
Eq.  7 is  difficult  because  the  integrand  is  oscillatory.  Further- 
more since  f^(x)  will  be  highly  skewed  to  the  right  of  its  maxi- 
mum, it  is  imperative  that  a quadrative  formula  be  employed  which 
yields  uniform  accuracy  for  large  values  of  x.  We  refer  the 
interested  reader  to  Appendix  A for  details  of  this  aspect  of 
the  problem. 


7 


On  the  principal  usee  of  the  characteristic  function  la  as 
a moment  generating  function  since 


mk  = <x 


k i 3k4> 

k>  = (-i)k  x 


at' 


t = o 


( 18; 


Unfortunately  Eq.  17  is  too  complicated  for  this  purpose,  but  it 
is  of  no  great  importance  since  the  moments  of  the  lognormal  PDF 
are  well  known.  The  moments  about  the  origin  are 


mk  * ek^  ek2o2/2 

and  Increase  very  rapidly  with  k.  The  central  moments  are 

nk  - EL(x-mi )k3  , 

and  Wartmann  has  shown  that 


(U;)k/2  k 


(19) 


(20) 


“k  = T , k/2~  I (-1)J  (?)  u(k-J)(k-j-l)/2 

(0)-l)n/<:  j =0  J 


where  u>  = exp  o2  and 


(21) 


We  will  utilize  these  formulae  shortly. 

The  coefficient  of  skewness  y,  and  the  coefficient  of 
excess  yz  are: 


8 


> 


(23) 


Y 


1 


^ 3 _ co  ^ (o)-l) 2 (g)+2)e 
lif2  a)^(a)-l)?/2e3^ 


(co-1)1'5  (a>+2) 


Y = ^ **  ^ _ to 2 ( to—  1) 2 (co‘*  + 2co3  + 3to 2 — 3 ) g **  ^ 

2 ” M * to 2 ( co—  1 ) 2e^ 

= a)4  + 2<o 3 + 3<o2  - 6 . (24) 

Neither  y i or  Y2  depends  on  £.  Furthermore  > 0,  thereby  indicat- 
ing a long  tail  to  the  right  of  the  maximum  of  the  lognormal  PDF. 
Some  numerical  values  of  Yj  and  y2  are  given  in  Table  2. 
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3.  NONUNIQUENESS  OF  MOMENT  PROBLEM  FOR  LOGNORMAL  PDF 

The  determination  of  a PDF  via  its  moments  (provided  they 
exist)  is  a time  established  method.  In  other  words  does  a set 
of  moments  determine  a PDF?  The  known  answer  is  that  at  least 
one  PDF  is  determined,  but  there  may  be  others  (i.e.,  the  solution 
may  not  be  unique) . There  is  a tendency  among  practical  users 
of  the  moment  method  to  assume  that  if  the  moments  exist  then 
they  determine  a PDF  uniquely,  such  is  not  the  case.  In  point  of 
fact  the  lognormal  PDF  cannot  be  uniquely  defined  by  its  moments'. 

The  Carleman  criterion2 for  the  uniqueness  of  the  PDF  defined 
over  ( 0 ,°°)  via  a moment  sequence  is  that 


-1/k 


CD 


(25) 


in  other  words  the  series  of  central  moments  must  diverge.  It  is 
not  difficult  to  show  that  increases  rapidly  enough  for  the 
series  in  Eq.  25  to  converge,  thereby  violating  the  Carleman 
uniqueness  criterion. 


This  indeterminate  nature  of  the  lognormal  moment  problem 
manifests  itself  in  two  important  ways- 

First,  any  theory  predicting  the  probability  density  of  the 
irradiance  fluctuations  of  an  optical  wave  which  has  as  its 
primary  goal  the  calculations  of  the  moments  must  be  checked  to 
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3oo  if  the  Carleman  uniqueness  criterion  is  valid.  In  particular 
any  theory  which  predicts  a lognormal  probability  density  for  the 
irradiance  fluctuations  strictly  on  the  basis  of  the  irradiance 
fluctuation  moments  must  be  suspect.  This  lends  credence  to  the 
contention  of  Klyatskin3 that  there  is  no  limiting  form  for  the 
probability  density  in  the  region  of  strong  fluctuations,  although 
the  criticisms  of  DeWol f4 concerning  the  macroscales  of  turbulence 
are  well  taken. 

Second,  there  is  the  very  interesting  and  provocative  paper 
by  Mitchell5 on  the  "permanence  of  the  lognormal  distribution." 
Mitchell's  thesis  is  that  the  sum  of  independent  lognormally 
distributed  random  variables  is  approximately  lognormal  in  ap- 
parent contradiction  to  the  central  limit  theorem  which  predicts 
a Gaussian  distribution.  He  employs  an  analysis  which  depends 
upon  the  moments  of  the  lognormal  PDF.  Since  the  moments  of  the 
sum  are  the  sum  of  the  moments,  the  Carleman  criterion  is  violated. 
Actually  there  is  no  contradiction  and  one  can  easily  derive  the 
main  points  of  Mitchell's  paper  by  standard  arguments  as  we  will 
demonstrate  in  the  next  section. 

As  an  historical  note  we  mention  that  Heyde7 was  evidently 
the  first  person  to  publish  the  fact  that  the  lognormal  is  not 
determined  by  its  moments.  His  proof  is  rather  involved  and 
does  not  employ  Eq.  25. 


11 


4. 


QUASIPERMANENCE  OF  THE  LOGNORMAL  PDF 


Let  us  now  consider  the  problem  of  determining  the  correction 
terms  to  the  central  limit  theorem  so  that  we  can  establish  the 
rate  of  convergence  to  the  Gaussian  PDF.  Fortunately,  Cramer7 has 
made  an  especially  detailed  study  of  this  problem  which  seems  to 
have  escaped  the  notice  of  many  workers  in  applied  probability 
theory,  yet  is  of  direct  importance  insofar  as  our  problem  is 
concerned.  Cramer's  analysis  is  fairly  general  and  we  will 
specialize  his  approach  to  the  case  where  all  the  random  vari- 
ables possess  exactly  the  same  PDF.  In  the  Cramer  approach,  the 
PDF  of  the  sum  of  the  independent  random  variables  is  given  as 

the  product  of  the  limiting  Gaussian  PDF  and  an  asymptotic  ex- 

-k 

pansion  in  powers  of  N . 

Let  us  return  to  Eq.  1 and  note  that  since  the  xk  are  inde- 
pendent random  variables,  then  the  mean  and  variance  of  the  sum 
x is 


<x>  = N£  , var  x = No2  . 


(26) 


It  is  convenient  to  rescale  the  random  variable  x by  defining 
new  random  variables  z^  and  z 


z 


k 


vi 

o 


z = 

a/N 


(27) 


so  that 
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<z>  = 0 


var  z = 1 . 


(28) 


The  characteristic  function  of  z is 


but 


E[eitz]  = 4 (t)  = [0  (t)]N 

z zR 


(t)  = 
zk 


°°  lZ|t 

e f(zk)dzk 


^ (it)n<(xk-£)n> 


n=0 


Tn  . ,,n/2 
o n!N 


hence , 


♦z(t>  = 


1 " 211  ' 


it  3m 


3 

r + 


t4p, 


N 


6o3N*  24o4N2 


+ 0(ir5/2) 


Taking;  the  logarithm  of  both  sides  yields 


-2  it3u 

log*  (t)  = p r 

z 2 6o  3lr 


h*(— ) 

* \24o4N/ 


t4  + 0(N"5/2) 


where  we  have  employed 


log(l  + s)  = s - ^ 


(s  < 1)  . 


Consequently , 
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(29) 


(30) 


(3D 


(32) 


(33) 


i. 


■ — - 


- — -- 


— ^ — — 


4>z(t)  = e 


= .-t2/2 


1 - 


iYlt: 

6n* 


Y2fc 


Yft1*  3 

+ 0(N"'f) 


"2^N  " “72N 


Note  that 


♦z(°)  = 1 , 


<Mfc>  = 4>«(-t) 


so  it  is  truly  a characteristic  function, 


(3^0 


(35) 


The  Fourier  transform  of  <f>2(t)  is  f (z) 


- r°o 

Vz)  = 27  *.(t)  e-lztdt  . 

— 00  " 


(36) 


Upon  substituting  Eq.  34  into  Eq.  36 


we  can  show  that 


r „(z)  - -J-  e-z!/2 
" /5¥ 


1 - 


h3(z)  + mt  h4(z)  + My)  + o(n"^) 


(37) 


where  the  hk(z)  are  the  hermite  polynomials  defined  in  Eq.  15  and 
Yl,Y2  are  the  coefficients  of  skewness  and  excess  respectively, 
as  defined  in  Eqs . 23,  24.  An  equivalent  expression  is 


fN(Z)  - *(2)  - ^ ^ + *(«)(2)  + 


0(N“*) 


(38) 


where 


<Kz) 


i-  e’z2/2 

2u 


(39) 


( n ) 

and  <f>  (z)  denotes  the  nth  derivative  with  respect  to  z.  This 

formula  is  very  useful  because  <f(z)  and  its  derivatives  are  ex- 
tensive tabulated.  Equations  37  and  38  are  special  cases  of 
Cramer’s  general  formula. 

The  leading  term  of  Eq.  37  is  the  dominant  term  as  N 
approaches  infinity  and  is  the  Gaussian  PDF  of  the  central  Unit 
theorem.  However,  for  large  finite  N,  the  other  terms  exert 
influence.  The  first  correction  term  proportional  to  N“^  contains 
the  skewness  coefficient  Yj  and  is  the  dominant  correction  term. 
Most  investigations,  such  as  for  sums  of  independent  random  sine 
waves,  involve  PDF's  that  are  symmetric  so  that  Yj  = 0.  In  that 
case,  the  degree  of  approximation  is  governed  by  the  coefficient 
of  excess  y2  which  decays  as  N"1 . Thus  sums  of  independent 
random  variables  having  nonzero  mean  converge  very  slowly  (i.e., 
as  N ) to  the  Gaussian  PDF.  This  is  the  situation  encountered 
when  the  random  variables  are  lognormally  distributed.  Further- 
more, when  y j r 0,  then  fN(z)  is  not  symmetric  since  h3(z)  = 
-h3(-z). 


In  order  to  calculate  the  skewness  coefficient  y1 (z)  for 
f N ( z ) , we  can  easily  show  that  the  third  moment  of  z is 
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(40) 


Consequently,  the  skewness  coefficient  for  fN(z)  is 
Vz>  - - ♦ O(N^) 


where  we  nave  employed  Eq.  23.  The  skewness  coefficient  for  the 
sum  decreases  only  as  and  the  tail  to  the  recultant  PDF  is 
very  noticeable  for  large  H.  The  reason  for  the  negative  sign 

lies  in  the  definition  of  the  normalized  variable  z since  z can 
only  be  positive  when  x > N£. 


Thus  fN(z)  is  a PDF  with  a long  tail  and  it  looks  as  if  it 
were  lognormally  distributed.  It  is  in  this  sense  that  the  leg- 

normal  PDF  is  "quasipermanent"  since  the  decay  to  the  limiting 
Gaussian  PDF  is  so  slow. 


APPENDIX  A 


The  evaluation  of  Eq.  7 is  approached  by  noting  that  <J>(t) 


essentially  vanishes  for  very  large  t;  in  Pig.  1 the  maximum 
value  of  t for  which  <J>(t)  is  nonzero  is  t = 20.  Thus  we  can 
write 


f(x)  = 


2tt 


T 


-T 


♦ <t)  e"ixtdt 


(A.l) 


where  T is  to  be  determined  according  to  the  specific  situation. 


For  large  values  of  |x|,  the  graph  of  the  integrand  consists 
of  positive  and  negative  areas  of  nearly  equal  size.  The  ad- 
dition of  these  two  areas  results  in  substantial  cancellations 
with  loss  of  accuracy.  In  addition,  a very  fine  mesh  spacing 
of  quadrature  points  is  necessary;  since  for  large  |x|,  the 
integrand  oscillates  very  rapidly.  Unfortunately  for  us,  it 
is  precisely  in  this  region  of  large  x that  reasonably  accurate 
values  of  the  integral  are  required.  The  method  developed  by 
the  author  is  based  upon  the  fact  that  f(x)  is  basically  a 
(finite)  Fourier  integral.  Now  in  the  advanced  mathematical 
literature  on  numerical  analysis,  there  is  a theorem  which 
roughly  states  that  the  trapezoidal  rule  of  numerical  integra- 
tion is  in  the  limit  of  very  email  meeh  size  the  way  to 
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evaluate  Fourier  integrals.  Unfortunately  the  mesh  size  must 
be  vanishingly  small  which  means  that  the  number  of  quadrature 

points  is  extremely  large  (tending  to  infinity).  This  fact 
mitigates  against  the  literal  use  of  the  trapezoidal  rule,  al- 
though the  Fast  Fourier  Transform  (FFT)  method  is  essentially 
a naive  version.  Nevertheless,  some  modification  of  the 
trapezoidal  rule  to  finite  mesh  sizes  consistent  with  o'  r ac- 
curacy requirements  is  possible.  Furthermore  we  really  need 
algorithm  which  retrains  uniform  accuracy  when  |x|  is  large. 

The  ordinary  trapezoidal  rule  when  applied  to  Eq.7  yields 

f(x)  = 2iM  V(m«t)e-lmx'5t  (A. 2) 

where  M6t  = T and  (2M  + 1)  ■ number  of  quadrature  points. 

The  weight  factors  are 

H = 1 , m i-  t m 

in 

■ 2 * m = ± M (A. 3) 

This  well  known  algorithm  is  based  upon  the  hypothesis  that 
the  entire  integrand  varies  linearly  over  the  mesh  interval 
6t  so  that  it  cannot  be  usefully  employed  unless  |x6t|«|. 

For  when  |x6t|  does  not  satisfy  this  stringent  requirement,  the 
integrand  varies  nonlinearly  over  6t. 


However  let  us  still  retain  the  trapezoidal  rule  but  modify 
it  by  requiring  that  only  <j)(t)  vary  linearly  over  the  interval 
fit.  In  this  manner  we  can  separate  out  the  strong  influence  of 
large  |x|  on  the  oscillations  of  the  entire  integrand.  In  a 
separate  paper  (to  be  published)  the  author  has  shown  that  un- 
der this  restriction  the  weight  factors  are  no  longer  constants 
but  depend  on  x: 


A 


H - 
m 


|sin  t;-  xfit 
xfit 


= 

1 - cos  xfit 

+ i 

xfit  - sin  xfit 

(xfit)2 

(xfit)2 

as 

-1 

1 - cos  xfit 

- i 

xfit  - sin  xfit 

1 (xfit)2 

(xfit)2 

m / ±M 


m = +M 


m = -M 

(A.H) 


Note  that  these  weight  coefficients  tend  to  those  of  the 
regular  trapezoidal  rule  as  (xfit)+0.  The  fact  that  only 
4>(t)  has  to  vary  linearly  over  fit  means  that  we  can  take  the 


mesh  size  to  be  relatively  coarse  in  many  cases  of  practical 
interest . 
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TABLE  1.  FIRST  TWELVE  a„  COEFFICIENTS 

n 

a0(s)  = 1 
(s)  = 0 
a2(s)  = s 
a3(s)  - s 
a4 (s)  = 3s2  + s 
a$  (s)  = 10s2  + s 
ag(s)  = 15s3  + 25s2  + s 
a?(s)  = 105s3  + 56s2  + s 
a8(s)  = 105s14  + 490s3  + 119s2  + s 
a9(s)  = 1260s4  + 1918s3  + 246s2  + s 
a10(s)  = 945s5  + 9450s4  + 6825s3  + 501s2  + s 
atl(s)  = 17325s5  + 56980s4  + 22935s3  + 1012s2  + s 
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TABLE  2.  NUMERICAL  VALUES  OF  y,  AND  y2  AS  A FUNCTION  OF 


0 

Yi_ 

li 

.1 

.30 

.16 

.2 

.61 

.68 

.3 

.95 

1.64 

1.32 

3.26 

.5 

1.75 

5.90 

.6 

2.26 

10.3 

.7 

2.89 

17.8 

.8 

3.69 

31.4 

.9 

4.75 

57.4 

1.0 

6.18 

110.9 
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w<t> 

Figure  1.  Characteristic  function  *x(x)  of  lognormal  probability 

density  function  with  o2  * 0.25i  real  part  of 

*x(t),  imaginary  part  of  $x(t). 
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CHAPTER  II 


REALIZABILITY  CONDITIONS  ON  THE  COVARIANCE  OF 
THE  WAVEFRONT  ABERRATION  FUNCTION 


f 


The  study  of  the  influence  of  random  wavefront  aberration 
functions  on  the  imaging  characteristics  of  optical  systems  is  a 
subject  of  current  interest.  The  primary  quantity  that  enters 
into  the  analysis  is  the  covariance  function  of  the  random  wave- 
front  aberration  function.  The  purpose  of  this  note  is  to  show 
that  the  class  of  admissable  covariance  functions  is  restricted 
to  chose  that  represent  an  isotropic  random  surface. 

The  complex  diffracted  amplitude  a(x,y)  in  the  specified 
receiving  plane  is 

a(x,y)  = jj  A„(p,q)  eikW(p,q)  ei(xp+yq)  dp(Jq  (1) 

aperture 

where  Afl(p,q)  is  the  amplitude  distribution  over  the  exit  pupil 
and  W(p,q)  is  the  wavefront  aberration  function  (Hamilton's 
mixed  characteristic) . Both  Afl  and  W are  real  functions  of  p,q. 

When  the  wavefront  W is  spatially  random  then  it  consists 
of  two  terms 

W(p,q)  = WD(p,q)  + Wg(p,q)  (2) 

where  WD  is  the  usual  deterministic  part  due  to  the  aberrations 
of  the  optical  system  and  Wg  is  the  stochastic  part  due  to  atmos- 
pheric turbulence.  Wg  is  taken  to  be  a spatially  stationary 
(homogeneous),  Gaussian  random  process  having  a zero  mean  and  a 
specified  covariance. 
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(3) 


E[Ws(p,q)]  = 0 


E 


= o2r(al-a2 ,6^6^  . 


(M 


Only  the  difference  of  the  independent  variables  enters  because 

W is  homogeneous.  In  the  vast  majority  of  cases,  W is  also 
s s 

isotropic  so  that  the  correlation  function  r becomes  a function 
only  of  p = [(pj-p2)2  + ( q 2 -q2 ) 2 ]^ . 


Since  W (p,q)  is  a surface  (it  is  the  surface  of  constant 
s 

phase) , then  its  slope  and  curvature  must  exist  and  these  require- 
ments place  a restriction  upon  admissable  forms  of  its  correla- 
tion function.  In  the  context  of  the  mean  square  stochastic 
calculus,  the  requirements  translate  into  conditions  on  the 
derivatives  of  the  isotropic  correlation  function.  Longuet- 
Higgins  (1)  has  shown  that  isotropic,  Gaussian  random  surfaces 
can  be  described  by  three  invariants  mfl,  m2 , and  m4  which  can  be 
expressed  in  terms  of  the  isotropic  correlation  function 


m 


d2nr( p) 


2n 


dp 


2n 


n = 0,1,2 


p = 0 


(5) 


Consequently,  in  order  that  the  three  parameters  exist,  the  iso- 
tropic correlation  function  of  W must  be  smooth  at  the  origin 


• 


P = 0 in  the  sense  that  it  and  its  second  and  fourth  derivatives 
exist . 


The  Gaussian  correlation  function 


employed  by  Barakat  and  Blackman  (2,3,4)  obviously  satisfies  the 
differentiability  conditions  since 


r"(0)  = (2a)2  ; r'"'(0)  = (2a)4  . (7) 

Not  so  with  the  Kolmorogov  five-thirds  spectrum  employed  by 
Roddier  and  Roddier  (5) 

5/ 

r(p)  - 1-  p . (8) 

Here 


r" ( p)  = - 
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T p 


r""(p)  = p"7/s 


(9) 


and  both  derivatives  become  infinite  (and  thus  do  not  exist)  at 
p = 0.  There  is  an  inherent  theoretical  contradiction  in  the 
Roddiers ' analysis  since  their  random  wavefront  does  not  allow 
slopes  and  curvatures  to  exist  though  they  proceed  to  obtain 
numerical  data.  The  reason  why  this  does  not  amount  to  a contra- 
diction in  practice  is  that  the  sampling  interval  is  finite;  the 
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effect  of  a finite  sampling  interval  being  to  filter  out  high 
frequency  (i.e.,  short  wavelength)  components  and  thus  alter  the 
behavior  of  the  covariance  function  at  the  origin. 


An  interesting  variant  of  these  conditions  holds  for  the 
slit  aperture,  here 


l < x ) - [+  A„(p)  elkW(P)  elxP  dp  . 

1 


Now  W(p)  being  one-dimensional  only  the  slope  condition  is  re- 
quired. This  means  that  the  permissible  forms  of  the  covariance 
function  of  W(p)  need  only  possess  finite  second  derivatives  at 

p = 0 . 


The  question  naturally  arises  as  to  what  happens  when 
AQ(p,q)  is  also  random.  Now  A0(p,q)  is  not  a surface  in  the 
sense  that  W(p,q)  is  and  there  are  no  additional  differentiability 
requirements  on  its  covariance  function. 


These  differentiability  conditions  become  much  more  compli- 
cated when  W(p,q)  is  nonisotropic  and  then  involve  nine  functions 
not  just  three  but  they  still  only  involve  the  equivalent  of  the 
fourth  moment. 


1 1 
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CHAPTER  III 

STATISTICS  OF  THE  TRANSFER  FUNCTION: 
AND  RANDOM  AMPLITUDE  EFFECTS 


TEMPORAL 
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1.  INTRODUCTION 


In  prior  work1 > 2 we  have  investigated  the  behavior  of  the 
transfer  function  of  an  optical  system  subject  to  random  wave- 
front  errors,  by  developing  explicit  formulae  for  the  moments  of 
the  transfer  function  random  process,  and  by  computing  the  low 
order  moments  in  addition  to  simulating  the  process  in  a digital 
computer.  Analysis  was  predominantly  restricted  to  errors  of 
phase,  i.e.,  effects  due  to  random  amplitude  errors  were  neglected. 
Also,  statistical  dependence  of  the  transfer  function  process  on 
time  (and  concurrently  on  integration  time)  was  not  part  of  the 
formulation.  As  the  formation  of  images  through  intervening  tur- 
bulence is  a present  and  practical  problem  at  hand,  considerations 
of  temporal  effects  and  those  of  random  amplitude  errors  should 
rightly  be  made,  it  is  the  inclusion  of  these  considerations, 
temporal  dependence  and  random  amplitudes,  within  the  framework 
of  our  prior  formalism,  to  which  this  chapter  is  addressed. 


2. 


TEMPORAL  EFFECTS 


Before  explicitly  treating  the  temporal  behavior  of  the 
transfer  function,  we  review  the  time-independent  formulation. 

To  aid  clarity,  only  one  spatial  (reduced)  coordinate  is  treated. 
Following  Barakat  and  Blackman2 , and  ignoring  here  random  ampli- 
tude errors, 


exp 

4-¥) 

(i) 

with  | a | < 2.  Equation  1 is  the  transfer  function  of  an  aberra- 
tion-free optical  system  with  a slit  aperture  operating  in  the 
Fraunhofer  receiving  plane.  The  aberration  function  thus  com- 
prises only  a stochastic  part,  taken  to  be  from  a real,  homogeneous 
Gaussian  random  process  having  zero  mean,  and  covariance  given  by 

- E[w(p + 1 °.)w(p  - 1 «,)] 

= °2r(  | otj -a2  j ) (2) 

where  0 < |r|  <1.  Here  and  throughout,  E[  ] is  the  expectation 
operator  with  respect  to  the  ensemble.  If  we  take  the  ensemble 
average  of  Eq.  1,  we  derive  the  well-known  result1*3 


T( 


«>  -i  ’ 
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E[T(  a)  ] = exp{-k2o2[l-r(a)] } - -L^ij  , |ct|  < 2 

= 0 |a|  > 2 . (3) 

Thus,  the  first  moment  behavior  of  T(a)  factors  into  a determin- 
istic part 

Td(c)  - 1 - (,,, 

and  a stochastic  component, 

Tr(a)  = exp{-k2o2[l-r(a) ] } . (5) 

Barakat1  has  shown  that  the  factorization  of  E[T(a) ] into  sto- 
chastic and  deterministic  parts  depends  on  the  homogeneity  of  the 
wavefront  errors,  and  occurs  for  the  homogeneous  case  because, 
after  commuting  the  operations  of  expectation  and  integration, 
the  integrand  does  not  depend  on  the  integration  variable  and 
can  be  taken  outside  the  integral.  It  is  further  shown  in  Ref.  1 
that  higher  order  averaged  products  of  T(a)  do  not  factor  into 
stochastic  and  deterministic  components.  Since  the  treatment 
presented  is  not  for  merely  T(a)  but  for  T(a,0),  i.e.,  the  two- 
dimensional  aperture  case,  we  might  expect  a similar  result  when 
one  spatial  and  one  temporal  dimension  are  included.  In  fact, 
the  basic  result  does  not  depend  on  the  number  of  dimensions 
over  which  the  transfer  function  integral  is  computed,  again 
provided  that  the  covariance  function  of  the  Gaussian  random 
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I 

f 


process  representing  the  wavefront  statistics  is  only  a function 
of  coordinate  difference. 


We  now  consider  the  time-dependent  transfer  function.  Let 
Tt(a)  be  defined  as 

x 1 - ¥ 

VO)  . ial)  eXP!kl|W(P  + l * *) 

- w(p  - | ,t)  J j dpdt  . (6] 

To  compute  E[TT(a)],  we  must  first  define  a new,  bivariate  wave- 
front  covariance  function.  Retaining  our  prior  notation, 

V lai  ~a2  I*  = E|^(p  + 2 ai  (p  “ 2 01 2 * 2 ) j 

= a2r( |at-a2 |,  |t,-t2|)  . (7] 

Taking  the  expectation  of  Eq.  6,  commuting  expectation  and  inte- 
gration, and  evaluating,  we  have 

x 1 - ¥ 

E[TT(a)]  c p?  [ [ exp{*k202[l-r(a,O)  J}dpdt 

J"  i(i  . 1“L) 

= exp{-k2o2[l-r(a,0)]}  (l  - ij-LJ  , |a|  £ 2 . (8; 


The  general,  formal  extension  for  two  spatial  dimensions  can 
be  performed  by  inspection: 

E[TT(a,S)]  = exp{-k20  2 [l-r(ot,8,0)]}  • TD(a,B)  (9) 

with  Td  the  deterministic  transfer  function.  Thus,  we  have  the 
remarkably  simple  result  that  for  the  stationary  zero  mean  Gaus- 
sian random  wavefront  considered,  the  expected  value  of  the  trans- 
fer function  does  not  depend  on  the  collection  (integration) 
time . 


Since  it  is  well-known  that  short  time  exposure  transfer 
functions  are  quite  different  from  transfer  functions  involving 
long  integration  times,  the  observed  difference  must  depend  on 
the  higher  moments  of  T^Ca). 


The  second  order  averaged  products  of  the  transfer  function 
at  one  frequency  ( | ot i | ® |a2|)  are  of  two  forms 


and 


36 


(11) 


These  expressions  are  evaluated  by  converting  the  product  of  two- 
fold Integrals  Into  a four-fold  integral  of  products,  commuting 
the  expectation  and  integration  operations,  and  using  the  expres- 
sion for  the  characteristic  function  of  sums  of  normally  distri- 
buted random  variables4.  After  performing  these  operations  we 
have 


E[|Tx(a)|2] 


exp{-2k2a2[l-r(a,0) ] } 

4t2 


x exp^+k2a2[2r(pi-p2 ,tj-t2)  - r(pl-p2+a,t|-t2) 

- r(pi-p2-a,ti-t2)])dp1dp2dt1dt2  (12) 


E[T2 ( a) ] = -i—  exp{-2k2a2[l-r(a,0) ] } 
1 4t2 


x exp{-k2o2[2r(pi-pz ,tj-t2)  - r(pi-p2+a,tx-t2) 

_ r(PJ-P2-a,tl-t2)]}dPidP2dtidt2  . 
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(13) 


As  in  the  time-independent  case,  the  two  expressions  are  identi- 
cal except  for  the  sign  change  in  the  exponential  of  the  integrand. 
Equations  12  and  13  can  be  used  to  compute  the  variances  of  the 
real  and  imaginary  parts  of  the  transfer  function;  the  covariance 
between  real  and  imaginary  transfer  function  components  is  seen 
to  vanish  because  Eq.  13  is  purely  real.  Thus,  dropping  explicit 
dependence  on  a and  t for  convenience,  and  using  Tn  and  T_  to 
denote  respectively  the  real  and  imaginary  part  of  T, 


E[|T|2]  = E[T2]  + E[T2]  (14) 

E[T2]  = E[T2]  - E[T2]  (15) 

var  Tr  = | E[ | T | 2 ] + | E[T2 ] - (E[T] ) 2 (16) 


var  Tj  = \ E[ |T|2]  - \ E[T2 ] . 


(17) 


For  higher  order  averaged  products,  we  again  note  the  formal 
equivalence  of  the  integrands  of  the  transfer  function  T (a)  and 
the  transfer  function  T(a,0)  for  0 fixed  at  zero: 


T(o,0)  = 


expjki  w|p  + ^ , q)-w|p-|,  qj  |dpdq  . (18) 

r»f.  nr’p  ' L * 


Aperture 
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Thus,  the  general  expression  for  averaged  products  of  the  trans- 


fer function1  derived  for  two  spatial  coordinates,  is  also  that 
for  the  case  of  one  spatial  and  one  temporal  coordinate,  apart 
from  changes  in  the  limits  of  integration,  and  in  normalization. 
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3.  AMPLITUDE  EFFECTS 

The  problem  of  the  effect  of  random  amplitude  errors  on  the 
transfer  function  has  historically  been  simplified  by  the  assump- 
tion that  the  logarithmic  amplitude  fluctuations  are  normally 
distributed.  Using  this  assumption,  in  addition  to  those  of 
stationarity,  and  independence  between  the  amplitude  and  phase 
fluctuations,  it  has  been  shown5*6  that  the  expected  value  of  the 
random  transfer  function  component  itself  factors  into  one  part 
due  to  phase  errors  and  a second  part  due  to  amplitude  errors. 

In  the  analysis  below,  we  also  consider  the  phase  errors  and  log 
amplitude  errors  to  be  from  zero  mean,  stationary  Gaussian  random 

processes,  but  we  do  not  require  that  these  processes  be  inde- 
pendent . 

Let 

i - ifl 

T(a)  = ^ -(i  - ifi)  A°(p + ^ A°(p " expMw(p + * 

- w(p  - f)  jdp  • (19) 

Setting  ip  = Zn  (Ag/A) , A chosen  to  normalize  T(a),  we  have 
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(20) 


T(o) 


To  compute  E[T( a) ] , we  may  use  the  general  characteristic  function 
formula,  evaluated  for  n = 4 terms.  Somewhat  simpler  bookkeeping 
arises  if  we  combine  terms  and  use  the  formula  for  n = 2 terms. 


With 


Zj  = kw(p  + - i^(p  + |) 

Zj  - kw(p  - |j  + i*(p  - |) 

a,  - 1 . a,  - -1  ; B,  = P + | . B,  - P - f 

and 

/ 2 2 

E[exp(ialzi  + ia2z2)]  = exp L ^ ? ? a£am  ^ ^ ^ ( 

I 1=1  m=l 

(21) 

there  are  four  (=  n2)  terms  in  the  sum.  These  are 
l = m = 1 

k2E[  W2  ( p + f)]  - E[^2  ( p + |)]  - 2iE[  kW(  p + f)  Mp  + f)  ] 


(22a) 
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A = m = 2 


k2E[W2(p  - |)]  - E|>2(p  - |)]  + 2iE[kW(p  - -|)  p - “)] 


(22b) 


and,  for  A = 1,  m = 2;  A = 2,  m = 1 

-{ k2E[W(  p - |)  W(p  + |)]  + E[^(  p - £)  u,(p  + “)] 

- lkE[W(p  - |)  i|,(p  + “)] 


+ ikE[W ( p + |)  ,p(p  - “)]}  . 


(22c) 


Simplifying, 
A = m = 1 


k2a2 

w 


- of  - 


21kCT,..0,,r  ,.(°) 


■ — » — i-  . j.  . ' 

ip  w ip  wip 


A = m = 2 


k2a2  - o2  + 2ika  o.r  ,(0) 
w ip  w ip  wipx  ' 


A = 1,  m = 2;  A = 2;  m = 1 


(23a) 


(23b) 


■{k2a2r 
w w 


(a)  + 


a?r 


<P 


(a)  +WV*[rw*(«)  - 


( 23c) 


Here,  the  variances  of  the  phase  and  log  amplitude  errors  are 
given  respectively  by  a*  and  a2.  The  two  autocorrelation  functions 
and  the  cross  correlation  function  are  also  subscripted  with  their 
pertinent  variables.  Combining  terms  as  dictated  by  Eq.  21  and 
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performing  the  integration  indicated  by  Eq.  20,  which  again  does 
not  depend  functionally  on  p,  we  have 

E[T(  a)  ] = exp{-k2a^[l-rw(a)]  + ajtl+r^a)  ] -kia^O^a) 

- r*w(-a)]1  * i1  - ■¥)  • ^ 

Equation  2 4 is  just  the  expression  usually  derived  under  the 

assumption  of  independence  between  w and  with  the  exception 

of  the  term  involving  the  cross  correlation  function  r ,(a). 

w\p 

Clearly,  if  rw^(a)  is  zero  for  all  values  of  ot,  this  term  vanishes. 
But  if  rw^  is  identically  zero  throughout,  then  w and  ip  are  un- 
correlated and,  by  virtue  of  their  being  from  stationary  zero 
mean  Gaussian  random  processes,  they  are  also  independent.  How- 
ever, we  note  it  is  sufficient  but  not  necessary  for  rw(j,(°0  to 
be  zero  for  the  cross  term  to  vanish;  we  require  only  that 
r^ia)  - for  all  a.  Thus,  within  the  stationary  zero 

mean  Gaussian  assumption,  the  expected  value  of  T(a)  factors  into 
a product  of  four  components:  one  standard  deterministic  part, 

one  part  due  to  phase  errors  acting  alone,  one  part  due  to  ampli- 
tude errors  acting  alone,  and  a cross  term  which  is  of  unit 
modulus  and  which  depends  on  the  odd  part  of  rw(^(a)  for  its 
existence . 
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